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We report Hartree-Fock (HF) based pseudopotentials suitable for plane-wave calculations. Unlike 
typical effective core potentials, the present pseudopotentials are finite at the origin and exhibit rapid 
convergence in a plane-wave basis; the optimized pseudopotential method [A. M. Rappe et al. , Phys. 
Rev. B 41 1227-30 (1990)] improves plane-wave convergence. Norm-conserving HF pseudopotentials 
are found to develop long-range non-Coulombic behavior which does not decay faster than 1/r, and 
is non-local. This behavior, which stems from the nonlocality of the exchange potential, is remedied 
using a recently developed self-consistent procedure [J. R. Trail and R. J. Needs, J. Chem. Phys. 
122, 014112 (2005)]. The resulting pseudopotentials slightly violate the norm conservation of the 
core charge. We calculated several atomic properties using these pseudopotentials, and the results 
are in good agreement with all-electron HF values. The dissociation energies, equilibrium bond 
lengths, and frequency of vibrations of several dimers obtained with these HF pseudopotentials and 
plane waves are also in good agreement with all-electron results. 

PACS numbers: 02.70.Ss, 71.15.-m, 31.25.-v 



I. INTRODUCTION 

In nearly all plane- wave density functional calculations 
[l|, the use of pseudopotentials is essential in order to 
eliminate the atomic core states and the strong potentials 
that bind them. This results in replacing the electron- 
nuclear Coulomb interaction by a weaker (often angular- 
momentum dependent) potential. The advantage is the 
reduced computational cost, due to reduction in the num- 
ber of electrons and in the plane- wave basis cutoff 0]. 

In the physics community, pseudopotentials are typi- 
cally constructed using density functional theory (DFT) 
due to its success in electronic structure calculations 
[Ij 0j [1] ■ Pseudopotentials based on methods other than 
DFT are less widely available. One of the key ingredi- 
ents of most of these pseudopotentials is their "softness," 
meaning that they provide rapid convergence in a plane- 
wave basis. On the other hand, in the chemistry com- 
munity, Hartree-Fock (HF) pseudopotentials or effective 
core potentials (ECPs) are mostly used [1, 0, H, H, fiot - 
Most of the available ECPs, which are typically expressed 
in a Gaussian basis, are of limited use in plane-wave cal- 
culations, mainly because of their singular behavior at 
the origin and slow convergence in a plane- wave basis. 

Recently, there has been a growing interest in devel- 
oping soft-core HF pseudopotentials [nl, [H, [H, [3 . HF 
pseudopotentials have been applied in diffusion Monte 
Carlo calculations, and the results suggest that they are 
better suited for correlated calculations than DFT-based 
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pseudopotentials [Uj. In addition, HF pseudopotentials 
could be useful in certain calculations where a nega- 
tively charged reference state is needed; HF tends to 
bind electrons more strongly than DFT. Also, hybrid 
exchange-correlation functionals, which include some 
proportion of HF exact exchange (notably B3LYP |15| ) 
have proved successful in electronic structure calcula- 
tions, and whether HF or B3LYP based pseudopotentials 
would be useful in these type of calculations remains to 
be explored. 

Trail and Needs [l3| have recently found that norm- 
conserving HF pseudopotentials constructed using stan- 
dard norm-conserving pseudopotential methods develop 
a non-decaying and non-local tail. These pseudopoten- 
tials would generally lead to erroneous results, and even 
infinite energies in solid state calculations. Previously, 
pseudopotentials based on exact exchange within the op- 
timized potential method (0PM) were also found to de- 
velop a similar, but local and eventually decaying long 
range structure [3, [13 , which led to erroneous re- 
sults in applications. 

Different schemes have been proposed to cure the long- 
range behavior in the HF and 0PM type pseudopoten- 
tials mentioned above. Originally with the 0PM method, 
different groups developed procedures for reducing or 
eliminating the unphysical long-range behavior, while re- 
taining desired eigenvalue matching and norm conserva- 
tion of the pseudopotential [l^, ■ This procedure was 
later implemented in a self-consistent manner by altering 
the pseudo-orbitals over all space, which led to a small 
violation of the norm conservation of the core charge. 
Nevertheless, the pseudopotentials proved transferable in 
both atomic and solid-state calculations [l^. Trail and 
Needs fixed the long-range unphysical behavior using a 
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similar self-consistent approach [13| that also slightly vi- 
olates norm conservation. Ovcharenko et al. [12l | have 
recently developed HF-based pseudopotentials which are 
finite at the origin and are expressed in terms of Gaussian 
basis, but are still relatively hard for typical plane-wave 
calculations, as we show in Sec. IV. 

In this paper, wc construct soft-core HF based 
pseudopotentials following the Rappc-Rabc-Kaxiras- 
Joannopoulos (RRKJ) method [3]. We use a self- 
consistent approach to restore the Coulombic tail in a 
manner similar to that used by Trail and Needs [l^, [l^ . 
Pseudopotentials are highly non-unique entities, and it 
is useful to have several alternative forms available with 
different properties. For example, one of the advantages 
of the RRKJ construction scheme is the softness of the 
resulting pseudopotentials. 

The rest of the paper is organized as follows. We 
first review the general construction scheme of the RRKJ 
pseudopotentials, and the procedure we used to "local- 
ize" the pseudopotentials (remove the long tail behavior). 
The crucial steps in the validation of the pseudopotentials 
are the study of their convergence and transferability, i. e. 
how well the pseudopotentials converge in a plane-wave 
basis set, and how well they perform in environments 
different than the reference configuration. These steps 
are explored in Sections III and IV. In Sec. Ill, we study 
plane- wave convergence of the atoms, and we investigate 
the transferability by looking at the ionization energy, 
electron affinity, and excitation energy of the first- and 
second-row elements. In Sec. IV, wc study several dimers 
at the HF level, and compare both the all-electron and 
pseudopotential results. Finally, in Sec. V, we conclude 
by summarizing the main points. 



VHmr)}]=Yl J '^^ 



n'l' 



, $„z(r')$:,,,(r') 



(3) 



$„'i'(r). 



(4) 

For convenience, equal and opposite self interaction terms 
are included in both the Hartree and exchange terms. 
Therefore, the Hartree potential Vff[{(/)(r)}] = Vff[/9(r)] 
is a functional of the spherically averaged total electronic 

density p(r) = X;„i I (/-n/WP- 

The second step in the construction of norm-conserving 

pseudopotentials typically begins with the generation of 
a smooth set of pseudo-orbitals 0„i (r) to replace the AE 
orbitals 0„i(r), such that 0„i(r) equals the all-electron 
orbital beyond some cutoff distance Tc- Different norm- 
conserving pseudopotentials mainly differ in the form of 
4'ni{T) for r < Tc, and the guiding principles used to 
construct (r). 

In the RRKJ method, the pseudo-orbitals in the core 
region are expressed as a linear combination of spherical 
Bessel functions ji{qkr) 



k=l 



(5) 



The Bessel wave- vectors qk are chosen such that, 
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II. PSEUDOPOTENTIAL CONSTRUCTION 

Wc follow the usual density functional theory approach 
for the construction of the norm-conserving pseudopoten- 
tials starting from an all-electron (AE) atomic calculation 
0. However, in this case, we solve the Hartree- Fock in- 
stead of the Kohn-Sham equations. Our HF solver is 
adapted from the code of Ref. [l^. 

The pseudopotential construction starts by choos- 
ing an electronic reference state, and then solving the 
Schrodinger equation of the system for the eigenstates 
$„i(r) and eigenvalues e„;. The wavefunction can be 
written as $„;(r) = <t>ni{r)Yim{0 , (f) / r where Yim{0,(j)) 
are the spherical harmonics. The orbitals 4'niij') satisfy 
the radial Schrodinger equation: 



Vuv[{cp{r))] ^niir) = e„i(f>„i{r), (1) 



where f = -d^/{2dr^) + + 1)/ {2r^), t>on(r) = ~Z/r is 
the bare nuclear potential, and Z is the atomic number. 
Vhf[{0('')}] is the HF potential defined such that, 

VkF[{0(r)}]0„i('-) = VHmrmni{r)+K'mr)mni{r) 

(2) 



where [qur) and 4>ni (^) ^^'^ ^^'^ derivatives of ji [qur) and 
4'ni{T) with respect to r, respectively. Njy is the number 
of Bessel coefficients, Cnik- 

The key ingredient of the RRKJ method is the opti- 
mization of the Bessel coefficients, Cnik, to minimize the 
residual kinetic energy defined as follows: 



AT„i({c„a-},gc) = - / rfr$„,(r)V2$„i(r) 

'''dqq2||.„,(q)|2. 



where $„i(q) is the Fourier transform of real space 
pseudo-orbital $„;(r) = (f>ni{r)Yi,n{0,(l))/r, and is the 
target wavevector above which the contribution to the 
residual kinetic energy is as small as possible. For a par- 
ticular 7 the planewave energy cutoff required to achieve 
total energy convergence of AT„i is given by the square 
of the maximum value of qc over all orbitals [J ■ 

The optimization of Ar„/({c„/fc}, gc) is constrained to 
pseudo-orbitals that arc continuous at r^, with continu- 
ous first and second derivatives [in fact, one of these local 
constraints would be already satisfied due to Eq. ©J , and 
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that satisfy norm conservation of the core charge: 



(8) 



This procedure defines all the Bessel coefficients in 
Eq. ([5]), and Nb must be at least 3 to allow the constraints 
to be satisfied. Any additional parameters will be opti- 
mized by minimizing the residual kinetic energy due to 
the non-linearity of Eq. ([7]). Throughout this study, we 
used iVb = 6. 

The screened pseudopotential V"^\{r) is defined as the 
potential that makes the desired pseudo-orbital 4'ni{i') an 
eigenstate of the one-electron Schrodinger equation, with 
the same eigenvalue e„; as the corresponding all-electron 
state: 



(9) 



This equation can now be inverted, thanks to the node- 
less character of (pni (r) , to obtain the screened pseudopo- 
tential: 



1) 



2r2 



0„i(r) . (10) 



Note that Kjcr ('') ^ continuous function because of the 
continuity requirements imposed on the pseudo-orbital 
4>ni{,r) and its first two derivatives. 

The last step in the construction of pseudopotcntials is 
to remove the screening effects of the valence electrons, 
and to obtain the ionic pseudopotential. This is done by 
subtracting the Hartree and exchange-potential contribu- 
tions of the valence electrons from the screened potential 
V;^^(r). That is. 



(i>niir) 
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where {(f>{r)}y includes only the valence orbitals. Note 
that the last term in Eq. pT|) . the exchange potential, is 
orbital dependent, and its explicit construction is feasible 
because the pseudo-orbitals are nodeless. 

Up to this point, the pseudopotential construction 
recipe mirrors exactly the procedure done in regular 
norm-conserving pseudopotcntials based on DPT, except 
that the original Hamiltonian is the HE and not the KS 
Hamiltonian. However, the resulting ionic potentials, 
T^"jj(r), have different behavior for large r, which in turn 
stems from the corresponding largc-r behavior of the HE 
and KS orbitals. ft was shown by Handeler et 
that the HE orbitals decay at infinity as 
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where a — (— 2eHo)^^^ is determined by the eigenvalue 
of the highest occupied orbitals ehOj and is an or- 
bital dependent quantity 2lJ. This behavior is to be 



TABLE I: A summary of the cutoff radii (in atomic units) 
for the pseudopotcntials for s, p, and d channels used in the 
atomic calculations. The last column shows g™^", the largest 
gc value over all the orbitals (in Ry^^'^). The cutoff energy for 
a planewave basis is approximately the square of g™'"'. 
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contrasted with the exponential decaying behavior of the 
KS orbitals, (t)ni{r) oc cxp[(— 2e„i)^/^r] for large r. 

In DET-bascd pseudopotcntials, the exchange poten- 
tial is replaced by the exchange-correlation functional of 
the electron density p(r), which decays exponentially to 
zero for large r. Thus, the ionic potential l^on(^) — 
— Zy/r for large r where Zy is the valence charge den- 
sity. 

The effect of the special long-range behavior of the 
HE orbitals on the ionic pseudopotential V^on(^) "^^^ be 
understood by writing Vj^j^ir) as. 



Cn(0 = -- + VH[p-p:„] 



<t>ni{r) 



for r > Tc 



(13) 



where the Schrodinger equation of Eq. (1) is used. Here 
p[r) is the total electronic charge density while /5„(r) is 
the pseudo-valence charge density. In the HE case, the 
exchange potential is a non-local functional of the HE 
orbitals as shown in Eq. In particular, the lead- 

ing behavior of the numerator for large r 
where (3 



max„; Pni , while the denominator will go as 
. Therefore, the exchange potential, as com- 



puted in Eq. (|13p . is not guaranteed to decay faster than 
l/r. In fact, the exchange potential can even grow with- 
out bound for large r when computed in this way. This 
is the source of the non-local long-range non-Coulombic 
tail in the HE pseudopotcntials [13l |. 

In this work we follow the Trail and Needs self- 
consistent approach to localize the ionic pseudopotential 
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FIG. 1; (Color online) The upper panel shows the Ne po- 
tential Vion constructed using the RRKJ method before the 
tail correction. The inset shows AFion = V;",*, + Z^/r which 
displays the non-physical tail of the s-potential channel for 
large r. The lower panel shows V|oc-ion after the tail correc- 
tion procedure. In the inset, we show the large-r behavior of 

AVloc-ion = + Z„/r. 
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FIG. 2: (Color online) The upper panel shows the all-electron 
(solid lines) and the localized pseudo-orbitals (dashed lines) 
of Ne constructed using the RRKJ method. The lower panel 
shows the difference in the pseudo-orbitals before and after 
the localization procedure (scaled by 10"'^). 



p^ . The localized ionic potential V|oc-ion('')' which is 
forced to behave as —Z^/r for large r, is defined as 



r < Hoc 



loc— ion 



(r) = 



. Z/r - VniPc)] + Vh(Pc) - Z/r r> n^c 

(14) 

where pdr) is the core charge density, l/VJ = Hoc/ 16 is 
the characteristic distance for the localized ^"c-ionl'^) 1° 



approach —Z^/r, and Hoc is the localization radius. For 
all pseudopotentials, the value of ri^c for each Z-channel 
is fixed to be the same as the corresponding rc- lni{r) is 
defined as Pni + Qnifir, Hoc) where Pni and Qni are param- 
eters which are chosen to localize the ionic pseudopoten- 
tial Vj'^l^{r) in a self-consistent manner, and /(r, Hoc) is a 
polynomial function whose explicit form depends on the 
pseudopotential construction method (see below). 

The self-consistent procedure is performed as follows. 
First, starting with an initial guess for the p„i and q„i pa- 
rameters, the HF equations with a new ionic potential, 
as defined in Eq. (fTi|) , are solved for a new set of orbitals 
and eigenvalues. The Pni and qni parameters are then 
adjusted to minimize the error in the eigenvalues and log- 
arithmic derivatives [2^. The procedure is then repeated 
until Pni and Qni parameters are found that give eigen- 
values and logarithmic derivatives that closely match the 
all-electron values to the desired tolerance. 

Finally, we define /(r, Hoc) used in 7ni(r). For the 
case of Troullier-Martins pseudopotentials, we chose the 
same form as was used in Ref. [l3|, namely f{r,r\oc) 

(1 - 2r2/3rf^^) for r < Hoc, and /(r,noc) 
otherwise. For RRKJ pseudopotentials, we investigated 
a few different forms and found the results to be insen- 
sitive to the choices. We will report our results using 
f{r, Hoc) = 1 -r/2rioc for r < Hoc, and /(r, Hoc) = noc/2 
for r > Hoc. 
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FIG. 3: Amount of the violation of the core norm- 
conservation for the s, p, and d pseudo-orbitals (in units of 
10^^ electrons), respectively. There is less than 0.001 elec- 
tron shift between the valence and core regions. Pseudopo- 
tentials are constructed using the RRKJ (open circles), and 
TM (closed circles) methods. 

To illustrate the pseudopotential localization proce- 
dure, we show in the upper panel of Fig. [1] the pseu- 
dopotentials Vj^l{r) of Ne as constructed using the RRKJ 
method with the parameters given in Table I. In the lower 
panel, we show the same Ne pseudopotentials after the 
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TABLE II: Ionization energy study of the first- and second-row elements. The first column shows the all-electron ionization 
energy, while the other columns show the deviation from the AE result. Columns 2 and 3 show the pseudopotential ionization 
energies as obtained using the original and the localized RRKJ pseudopotentials [Eqs. (|I3|I and (|14p . respectively]. Columns 
4 and 5 show the same values as obtained using the TM method. For comparison, we show also the values obtained by Trail 
and Needs The pseudopotential parameters are summarized in Table H] We show also the average abs. error for each 

pseudopotential with respect to the AE values. All energies are in Hartrees. 
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TABLE III: Similar to Table ITT] except we show the electron affinity of the first- and second-row elements. 
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self-consistent localization procedure. As can be seen, 
the differences between the upper and lower panels are 
rather small, and arc in the core as well as in the va- 
lence regions. Most of the changes are within the core 
region, however. The insets show the long-range behav- 
ior of Vioni''') + Zi,/r. Note, in particular, in the upper 
panel, the non-decaying tail of the s-potential, which can 



be shown [l3| to behave as a + b/r for large r where a 
and b are real numbers. 

In Fig. [51 we show in the upper panel the all-electron 
and pseudo-orbitals after the localization procedure. The 
lower panel shows the difference in the pseudo-orbitals 
before and after the localization procedure. Again, 
changes are found in both the core and valence regions. 
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FIG. 4: Similar to Fig. 3, except that we show the relative 
error in the logarithmic derivative at Tc from the all-electron 
values. The relative error is defined as |1 — Lioc/L\, where L 
and Lioc are the logarithmic derivative before and after the 
localization procedure is applied, respectively. The logarith- 
mic derivative is approximately conserved in the localization 
procedure. 



One obvious drawback of the localization procedure is 
that the core-norm will be different from the all-electron 
value, i.e. Eq. ([5]) will be violated slightly. This can be 
seen in Fig. [2] because the changes in the orbitals for 
r < Tc are all positive, and will integrate to a non-zero 
value. 

As mentioned before, the localization procedure con- 
serves the all-electron eigenvalues, and the logarithmic 
derivative. In practice, the logarithmic derivative is only 
conserved approximately. Figures [3] and |4] show a sys- 
tematic study of the violations of the norm-conservation 
and the change in the logarithmic derivative at due to 
the localization procedure for the first- and second-row 
elements. The study was done using both the RRKJ and 
TM pseudopotentials. From Fig. [31 we see that there 
is less than k, 0.001 electron shift between the core and 
the valence. The relative error between the logarithmic 
derivatives before and afte r ap plying the localization pro- 
cedure is shown in Fig.[3] [2^. As can be seen from both 
figures, both pseudopotential construction methods yield 
similar results. 



III. PSEUDOPOTENTIAL TESTS ON ATOMIC 
PROPERTIES 

In this section, we study the transferability of the 
pseudopotentials by looking at several atomic properties 
namely ionization energies, electron affinities, and exci- 
tation energies of the first- and second-row elements. All 
of these calculations are performed by solving the HF 



equations in real space with the same code [25| which is 
used to generate the pseudopotentials. 

We show results obtained using HF pseudopotentials 
constructed using both the RRKJ and TM methods. For 
the TM pseudopotentials, wc used the same code [25| 
which is used to generate the RRKJ pseudopotentials 
but following the TrouUicr-Martins construction scheme. 

The parameters for all pseudopotentials used in the 
study of the atomic properties are summarized in Ta- 
ble HI In order to aid in comparison to previous results, 
the reference configurations and the construction param- 
eters are the same as of Ref. [l^. In general, however, 
one has to choose these parameters according to the tar- 
geted applications. Moreover, unlike the TM method, 
the RRKJ method has two additional adjustable param- 
eters for each pseudo-orbital, namely Ni, and Qc- We set 
Nb = & for all pseudopotentials, and selected the qcS 
such that ATni k, 5 meV/electron for each orbital. As 
mentioned before, the energy cutoff required to achieve 
this level of energy convergence in the target calculations 
is approximately the square of the largest value used. 

In Table [TI] we present the ionization energies. For 
each element, we calculated the ionization energy using 
our RRKJ and TM pseudopotentials. These values are 
compared with the all-electron value (shown in the first 
column), and we only report the difference, AE'ion, from 
the all-electron energy. AE'ion is obtained using the orig- 
inal pseudopotential of Eq. (|13p with the unphysical tail 
behavior, while Ai^ion is calculated using the localized 
HF pseudopotential of Eq. pi)) . We show also for com- 
parison the same results as obtained using the HF TM 
pseudopotentials as reported in Ref. p^ . 

For both pseudopotential construction schemes, the 
agreement between Ai?ion and Ai^ion indicates that the 
self-consistent procedure that is used to localize the pseu- 
dopotential did not change significantly the original po- 
tential. It is important to stress that the effects of the 
non-Coulombic tail in atomic calculations are minimal, 
and this is why there is a good agreement between Ai?ion 
and Ai^ion- However, in any solid calculation with pe- 
riodic boundary conditions, the results obtained using 
pseudopotentials with the tail problem would be erro- 
neous. 

Also, the pseudopotential results in Table |TT] are in 
good agreement with the all-electron values for both the 
RRKJ and TM methods. On average the difference is less 
than 0.5 milli-Hartree (mHa) and the largest deviation is 
less than 1.5 mHa. This is a clear indication of the good 
quality of the pseudopotentials. Finally, the results ob- 
tained using the TM construction scheme are in excellent 
agreement with the equivalent results obtained by Trail 
and Needs [l3| . Any small differences could be attributed 
to the different grids used, or the slight differences in the 
self-consistent HF procedure. 

In Table IIIIl we show a similar study for the electron 
affinity of the first- and second-row elements. The elec- 
tron affinity is the difference in energy between the neu- 
tral atom and the negatively charged ion in their ground 
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TABLE IV: Excitation energy study of the first- and second-row elements. The first column shows the ground state configu- 
ration, and the second column shows the configuration of the excited states. The all-electron excitation energy is shown in the 
third column. The other columns show the deviation from the AE result, similar to Table HIl All energies are in Hartrees. 



Atom 
Ground state 


Excited state 


AE 


RRKJ 

APexc 


APexc 


TM 

APexc 


APexc 


TM" 

APcxc 


APcxc 


H ls'[-'S] 


2p^pP] 
3d' I'D] 


0.3750 
0.4444 


0.0000 
0.0000 


0.0000 
0.0000 


0.0000 
0.0000 


0.0000 
0.0000 


0.0000 
0.0000 


0.0000 
0.0000 


He Is^l^S] 


ls'2p'fP] 
ls'3d'(-'D] 


0.7302 
0.8061 


-0.0008 
-0.0009 


-0.0008 
-0.0008 


-0.0006 
-0.0007 


-0.0006 
-0.0007 


-0.0006 
-0.0006 


-0.0006 
-0.0006 


Li 2s^('S] 


2p'['P] 
3d'['D] 


0.0677 
0.1408 


0.0000 
0.0000 


0.0000 
0.0000 


0.0000 
0.0000 


0.0000 
0.0000 


0.0000 
0.0000 


0.0000 
0.0000 


Be 2s^[^S] 


2s'2p'[^P] 
2s'3d'[^D] 


0.0615 
0.2389 


0.0011 
-0.0002 


0.0011 
-0.0002 


0.0010 
-0.0002 


0.0010 
-0.0002 


0.0020 
-0.0001 


0.0020 
-0.0001 


B 2s'^2p^[^P] 


2s'2p'[^P] 
2s''3d'^D] 


0.0784 
0.2353 


0.0023 
-0.0006 


0.0022 
-0.0006 


0.0020 
-0.0005 


0.0019 
-0.0005 


0.0020 
-0.0005 


0.0019 
-0.0005 


C 2s^2p^[^P] 


2s'2p^^S] 
2s^2p^3d^pP] 


0.0894 
0.3402 


0.0060 
-0.0008 


0.0058 
-0.0008 


0.0054 
-0.0006 


0.0053 
-0.0005 


0.0054 
-0.0006 


0.0053 
-0.0005 


N 28^2/ [^S] 


2s'2p''[^P] 
2s^2p^3dM*^l 


0.4127 
0.4565 


0.0029 
-0.0009 


0.0030 
-0.0008 


0.0030 
-0.0006 


0.0030 
-0.0006 


0.0031 
-0.0005 


0.0031 
-0.0005 


O 2s^2p-^[^P] 


2s''2p'^3d'f D] 
2s'2p^^P] 


0.3809 
0.6255 


0.0001 
-0.0003 


0.0001 
-0.0002 


0.0001 
-0.0002 


0.0001 
-0.0001 


0.0000 
-0.0002 


0.0000 
-0.0001 


F 2s^2p^[^P] 


2s'2p^3d'[^F] 
2s'2p''[^S] 


0.5220 
0.8781 


-0.0004 
-0.0052 


-0.0003 
-0.0051 


-0.0003 
-0.0051 


-0.0002 
-0.0049 


-0.0003 
-0.0051 


-0.0002 
-0.0049 


Ne 2s^2p^[^S] 


2s22p^3dipP] 
2si2p'*3dip_D] 


0.6734 
1.7565 


-0.0008 
-0.0051 


-0.0007 
-0.0049 


-0.0005 
-0.0047 


-0.0004 
-0.0045 


-0.0006 
-0.0048 


-0.0004 
-0.0045 


Na 3s' fS] 


3d'^D] 


0.0725 
0.1263 


-0.0001 
-0.0002 


-0.0001 
-0.0002 


-0.0001 
-0.0001 


-0.0001 
-0.0001 


-0.0002 
-0.0001 


-0.0002 
-0.0001 


Mg 3s' ["-S] 


2s'2p'[^P] 
2s^3d4^_D] 


0.0679 
0.1843 


0.0008 
-0.0002 


0.0008 
-0.0002 


0.0007 
-0.0002 


0.0007 
-0.0002 


0.0007 
-0.0002 


0.0008 
-0.0002 


Al Ss^Sp^pP] 


3s'3p'[*P] 
3s''3d'^D] 


0.0858 
0.1441 


0.0008 
-0.0002 


0.0007 
-0.0002 


0.0007 
-0.0002 


0.0007 
-0.0002 


0.0008 
-0.0002 


0.0007 
-0.0002 


Si Ss^Sp^pP] 


3s'3p\^S] 
3s''3p'3d'[^F] 


0.0913 
0.2146 


0.0023 
-0.0001 


0.0022 
-0.0001 


0.0023 
-0.0001 


0.0022 
-0.0001 


0.0023 
-0.0001 


0.0022 
-0.0001 


P Ss^Sp^'pS'] 


3s'3p\^P] 
3s"3p"3dipP] 


0.3006 
0.3023 


-0.0004 
-0.0008 


-0.0003 
-0.0007 


-0.0003 
-0.0007 


-0.0003 
-0.0007 


-0.0004 
-0.0007 


-0.0003 
-0.0007 


S 3s^3p^pP] 


3s^3p^3dipP)] 
3si3p^pP] 


0.2672 
0.4260 


0.0014 
-0.0010 


0.0014 
-0.0010 


0.0015 
-0.0010 


0.0015 
-0.0010 


0.0013 
-0.0009 


0.0015 
-0.0009 


CI Ss^SpSpP] 


3s23p''3di[''P] 
3si3p^p5'] 


0.3733 
0.5653 


-0.0001 
-0.0020 


-0.0001 
-0.0020 


0.0001 
-0.0020 


0.0002 
-0.0021 


-0.0002 
-0.0018 


-0.0001 
-0.0018 


Ar 3s^3p'^['S] 


3s^3p^3d^ pP] 
3si3p'^3dipP] 


0.4824 
1.1597 


-0.0014 
-0.0028 


-0.0013 
-0.0027 


-0.0014 
-0.0026 


-0.0012 
-0.0024 


-0.0014 
-0.0026 


-0.0013 
-0.0024 


Average abs. error 
Maximum abs. error 




0.0012 
0.0060 


0.0011 
0.0058 


0.0011 
0.0054 


0.0011 
0.0053 


0.0011 
0.0054 


0.0011 
0.0053 



Ref. [T3|. 
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TABLE V: The parameters for the pseudopotentials used in 
the dimer calculations. The cutoff radii are in atomic units 
and g^"" is in Ry'/^ 



Atom Tcs 


Tc.p 


max 

Ted Qc 


N 0.91 


0.91 


0.91 10.6 


P 1.58 


1.58 


1.90 5.5 


CI 1.68 


1.68 


1.90 6.2 


TABLE VL Planewave basis cutoff ( 


in Ry) for nitrogen 


and phosphorus pseudopotentials for several pseudopoten- 


tials methods. The 


planewave basis 


cutoff is estimated 


using Eq. ^ for a 


residual kinetic 


energy of AT„i 


5 meV/electron, which 


can also be read from Fig. [S] 






Ec„t(Ry) 




N 


P 


RRKJ 


112 


30 


TM 


145 


50 


OAL 


210 


55 


SBKJC 


375 


90 



state configurations. The results agree to within 0.5 mHa 
with previously published HF electron affinities (26j . 

Again, these results show that the pseudopotentials 
have good transferability properties. Also, we note here 
that the HF pseudopotentials generally bind an extra 
electron more strongly than the DFT-based pseudopo- 
tentials. This is why the study of the electron affinity is 
feasible. 

We report the results of excitation energies in Table HVl 
This is the difference in energy between the ground state 
configuration shown in the first column and the two ex- 
cited state configurations shown in the third column. As 
before, the pseudopotential results are in good agreement 
with the all-electron values. The absolute average devi- 
ations from the all-electron values are consistent across 
the different pseudopotential schemes, and is less than 
1.2 mHa. The largest deviation is with carbon and is 
w5-6 mHa (the excitation energy between the '^P and 
states). 



IV. STUDY OF SEVERAL DIMERS 

In this section, we apply our HF pseudopotentials in a 
study of three dimers and compare them with all-electron 
calculations. The pseudopotential calculations are done 
using a HF planewave basis code [2^ , and the all-electron 
results are obtained using a Gaussian basis code [1^ . For 
comparison, we have also generated LDA (soj pseudopo- 
tentials using exactly the same parameters as those of the 
HF pseudopotential, and we compare them with their all- 
electron counterparts on an equal footing with the HF re- 
sults. In our all-electron and pseudopotential results we 
used the Perdew-Wang [s^l flavor of LDA. We did not 
include a non- linear core correction fNLCOfssl. [s^ in 
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FIG. 5: (Color online) Comparison of several ionic pseudopo- 
tentials for nitrogen and phosphorus. We show pseudopo- 
tentials obtained using the RRKJ, TM, OAL [H], and the 
SBKJC [13 methods. For the RRKJ and TM methods, the 
Tc values are the same as those in Table FVl We show only the 
s and p potentials. 



our LDA pseudopotentials for LDA/HF comparison pur- 
poses. The pseudopotential DFT calculations are carried 
out using ABINIT [3l| . The all-electron calculations are 
performed using the correlation-consistent cc-pV5Z basis 
set for N and cc-pV(5+d)Z for P and CI Hi]. The dif- 
ferences between the quadrupole and quintuple basis sets 
are negligible; e.g., in the binding energies the differences 
are less than 0.02 eV. We also verified these values using 
6-311++G(3df,3pd) basis sets. 

In Table fVl we show the cutoff radii we used to gener- 
ate the RRKJ pseudopotentials for this study. These val- 
ues, which are different from those of Table [H are chosen 
to accommodate the bond lengths of these dimers with- 
out core overlap when possible. In the case of P and CI, 
the d core cutoff radius is allowed to extend beyond half 
of the dimer bond length to make softer pseudopoten- 
tials. The atomic reference configurations are the same 
as before. 

Before presenting the molecular results, we comment 



9 




50 100 
E (Ry) 

cut ^ -I ' 

FIG. 6: (Color online) Same as Fig. (5] except that we show 
the residual kinetic energy convergence (logarithmic scale) for 
the potentials calculated using Eq. ((Tj. 



on the planewave convergence of the RRK J pseudopotcn- 
tials in comparison to other pseudopotentials or ECPs. 
We examined the pseudopotentials of nitrogen and phos- 
phorus as representatives of the first- and second-row el- 
ements. Figure E] shows the RRK J, TM, OAL [i3|, and 
27[ pseudopotentials in real space for the two 



SBKJC 

elements. We did not include a d-channcl in these plots 
because neither the SBKJC nitrogen pseudopotential nor 
the OAL potentials have a c?-channel. For the RRKJ and 
TM pseudopotentials, we used the cutoff radii as shown 
in Table El Note that the SBKJC ECP has a divergence 
at the origin, and is thus not suitable for planewave cal- 
culations. More generally, direct examination of pseu- 
dopotentials in real space gives little information about 
the planewave basis needed for converged eigenstates and 
eigenvalues. 

One way to monitor the size of the planewave basis 
needed is to look at the residual kinetic energy conver- 
gence Q of the pseudowavefunctions, which we show in 
Fig. [SI The residual kinetic energy is calculated as shown 
in Eq. ([7]) by varying (we plot it against which is 
the cutoff energy, for convenience). The planewave basis 



TABLE VII: Dissociation energy in eV, bond length in 
Bohr, and harmonic frequency bJe in cm~^ for several dimers 
as calculated using density functional LDA and HE theories. 
We compare both the all-electron and the pseudopotential 
results. For each method, we generated the RRKJ pseudopo- 
tentials using the same level theory i.e. LDA and HF, respec- 
tively. 



LDA 



HF 







AE 


PSP 


AE 


PSP 






11.61 


10.96 


5.10 


4.92 




R. 


2.07 


2.06 


2.01 


2.02 






2388 


2380 


2720 


2722 


P2 




6.25 


6.02 


1.70 


1.70 




Re 


3.57 


3.57 


3.50 


3.50 






794 


788 


912 


905 


CI2 


De 


3.62 


3.55 


0.84 


0.84 




Re 


3.74 


3.74 


3.73 


3.72 




LOe 


561 


559 


610 


613 



cutoff for each pseudopotential is determined by the po- 
tential channel requiring the largest cutoff energy. The 
estimated planewave basis cutoffs based on Fig. [5] are 
summarized in Tabic IVII As can be seen, RRKJ pseu- 
dopotentials give the smallest planewave basis cutoffs. 

We summarize the results for the spectroscopic proper- 
ties of the dimers in Table EIIl The LDA and HF results 
for the equilibrium bond length and the harmonic vibra- 
tional frequencies of the dimers are shown to be well re- 
produced by both types of pseudopotentials. In the LDA 
dissociation energies there is a large discrepancy between 
all-electron and pseudopotential results for N2, and to a 
lesser extent in P2 and even much less in CI2. The HF 
pseudopotential dissociation energy of N2 shows a signif- 
icant deviation from the all-electron HF result, whereas 
for P2 and CI2, the results are the same to within 0.01 eV. 
A similar difference between all-electron and pseudopo- 
tential results was also seen in ReL [l^ using the 0PM 
method. 

Most of errors in the LDA results arc due to linear 
descreening of the pseudopotential, which is typically re- 
paired by including a NLCC in the target calculation. 
Furthermore, because all three dimers are closed shell, 
the linear descreening effects are largest in the isolated 
atom calculations and therefore mostly affect the disso- 
ciation energies. As shown by Porezag et al. [s^], these 
errors arc larger for first row species compared to heav- 
ier atoms and increase with more unpaired electrons, ex- 
plaining why N2 has the largest error (first row, and 3 un- 
paired spins) and why CI2 (second row, 1 unpaired spin) 
has the smallest. Our LDA all-electron and pseudopo- 
tential energies are in excellent agreement with previous 
results [l^ without a NLCC. This study also showed that 
the all-electron results are in good agreement with the 
pseudopotentials values after including a NLCC, indicat- 
ing that this is the dominant error for LDA pscudopo- 
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tentials. 

At the HF level there is also an error in the exchange 
contribution to the pseudopotential since the terms aris- 
ing from core- valence interactions in the atomic reference 
state are frozen into the core and only valence- valence ex- 
change terms are subtracted during the descreening step. 
As shown by our results and others [bsI . [36j , this type of 
descreening error is much smaller in magnitude in HF 
pseudopotcntials compared to LDA. This would suggest 
that HF derived pseudopotcntials are advantageous over 
those from LDA, at least in better descreening of the 
core- valence contributions. 



V. SUMMARY AND CONCLUSION 

Constructing HF pseudopotcntials suitable for 
planewave calculations is a non-trivial task considering 
the extreme nonlocality of the exchange potential. Pseu- 
dopotcntials constructed using a typical norm-conserving 
procedure develop a non-local and a non-Coulombic 
tail for large r. This would generally lead to erroneous 
results, especially in solid calculations. Several schemes 
have been introduced to cure the long tail behavior 
which lead to a small violation of the conservation of the 
core charge [ll[ll[ll. 

In this study, we present soft-core HF pseudopotcn- 
tials constructed using the RRKJ procedure which opti- 
mizes the potentials to yield rapid planewave basis cut- 



off convergence. The long tail behavior is fixed using 
a self-consistent procedure following that of Trail and 
Needs [l^l, which leads to a negligibly small violation of 
norm-conservation. These pseudopotcntials are applied 
in HF calculations of several atomic properties yielding 
results in good agreement with the all-electron values. 
We also apply them in the study of the dissociation en- 
ergies, equilibrium bond lengths, and frequency of vibra- 
tions of several dimers, using a HF planewave code [28| . 
The all-electron and pseudopotential results are in agree- 
ment with each other, and the values are consistent with 
a similar comparison done using LDA pseudopotential 
and LDA all-clcctron calculations. 

Generation of HF based pseudopotcntials has been re- 
leased in version 3.0 of the GPL package OPIUM [Hi 
and is available for download. 
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